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Abstract 

We calculate the contribution to the elliptic flow observable f 2 from two-particle cor- 
relations in minijet production in ultrarelativistic heavy ion collisions. We use a minijet 
production cross section derived in a model inspired by saturation approach to high en- 
ergy scattering. Resulting differential elliptic flow V2{pt) is an increasing function of pT 
for transverse momenta below the saturation scale Qs- At higher transverse momenta 
{pt > Qs) differential flow stops growing and becomes approximately constant, repro- 
ducing the elliptic flow saturation data reported by STAR. The centrality dependence of 
the minijet contribution to V2 is also in good agreement with the data. 
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1 Introduction 



Ultrarelativistic heavy ion collisions provide us interesting new information on QCD under 
extreme conditions. Different stages of the collisions probe QCD in different regimes: early 
times are likely to be characterized by strong saturated gluonic fields while creation of 
quark-gluon plasma (QGP) is an attractive possibility in the later stages of the collisions 0. 
To be able to understand the experimental data generated in the collisions one has to learn to 
disentangle between the contributions of these two (in principle) different physical processes to 
heavy ion observables. 

The strong quasi-classical gluon fields are produced in the early stages of the collisions due 
to high transverse densities of color charges in the wave functions of the Lorentz-contracted 
nuclei p|, ^, |[ ||. The typical color charge density is characterized by the saturation scale 



Q^, which grows with nuclear atomic number and with energy of the collision [1^, [11 

While being coherent over large longitudinal distances these quasi-classical fields have a rather 
short coherence length in the transverse direction, of the order of the inverse saturation scale 
1/Qs- Saturation-inspired models have been quite successful in describing the emerging RHIC 
data on rapidity and centrality dependence of the total charged multiplicity of the produced 
particles [0 0. 

In the later stages of the collisions gluons and quarks generated by the quasi-classical fields 
are likely to reach kinetic, and, possibly, chemical thermal equilibrium, producing the quark- 
gluon plasma (QGP) 0, |14|, |l^, |l^. To be able to study this new state of matter one should be 
able to distinguish the contributions to physical observables of the collective phenomena due 
to QGP from the effects of strong gluonic fields of the early stages. Ideally one should try to 
construct observables which would be sensitive only to one type of physics while being inde- 
pendent of the other. An example of such observables are the long range rapidity correlations 
predicted in |]T^, which are almost entirely due to the dynamics of initial conditions. 

In this paper we would like to study the contribution of the quasi-classical gluonic fields to 
the elliptic flow observable V2, defined as the second Fourier moment of the azimuthal momen- 
tum distribution of the produced particles [18|. Elliptic flow reflects anisotropy of the transverse 
momentum distribution of the produced particles with respect to reaction plane. There has 
been a large amount of elliptic flow data produced in the heavy ion collisions at SPS ||T9|, EOf and 



RHIC |2l|, |2^. While the majority of these data are in good agreement with hydrodynamic 
simulations p4| , p5[1 , the emerging new STAR data on differential elliptic flow V2{pt) at high pt 
seems to deviate from hydrodynamic predictions |26|. Instead of continuing to increase with pt 



the flow variable V2{pt) saturates to approximately a constant above pt = 1.5 GeV E^. The 
data goes up to pt ~ 4.5 GeV [|2^. It is natural to expect that at such high momenta the hard 



(perturbative) physics should be responsible for the underlying strong interactions dynamic^. 
This led the authors of to propose a model of interplay of soft and hard interactions in V2'- 
while the low momentum part of f 2 was still described by hydrodynamic calculations the high 
Pt part was described by medium-induced radiative energy loss of partons ||2^ which would 
generate azimuthal anisotropy in momentum space due to coordinate space anisotropy of the 
overlap region. The resulting flow observable would saturate at some relatively large px and 
then would decrease with increasing pt- 

*Very recently an attempt was made in p9|] to explain the high-pT elliptic flow using pure classical gluon 
fields of the nuclei 0] yielding the elliptic flow which is too small to explain the data. 
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Here we are going to propose a model of a non-flow particle production mechanism gener- 
ating an elliptic flow observable f 2 which would go to a constant at large pt and would stay 
approximately constant as pt increases. Let us picture particle production at the early stages 
of a nuclear collision at high energy. In the first approximation high-p^ particles are produced 
independent of each other. We can illustrate this in the framework of the saturation approach 
to nuclear collisions: there the gluon fields are coherent over very short transverse distances 
of the order of l/Qs, with Qg ^ ^qcd- The high energy wave function of a large nucleus 
can be viewed in the transverse plane as consisting of many independent gluon fields each of 
them occupying a transverse area l/Q^. (In the longitudinal direction each classical field is 
of course coherent over the whole nucleus p| and this way the nucleus can be considered as 
"sliced" into narrow pipes with diameter l/Qg.) Therefore a collision with the overlap trans- 
verse area of the two nuclei S± could be pictured as at least S±Ql independent (sub) collisions. 
For sufficiently central heavy ion collisions S±Ql ~ Npart 3> 1, where Npart is the number 
of nucleons participating in the collision. At the leading order in this large parameter S±Ql 
the two-particle multiplicity distribution factorizes into a product of two single-particle mul- 
tiplicity distributions. Correlation between a pair of particles in the perturbative production 
mechanism happens when both particles are produced in the same subcollision, i.e., at the same 
impact parameter. Therefore these correlations appear at the subleading in S^Ql order and 
are suppressed by a power of S±Ql. We are going to argue below that this does not prevent 
them from significantly contributing to elliptic flow observable. 



The standard definition of differential elliptic flow is [pO| , |3^, |35 



V2{pt) = ( e'^^f^T-'^-^ ) = { cos2(0p, - ) (1) 



where (pp^ is the azimuthal angle of the produced particle with the value of transverse momen- 
tum Pt, $_r is the azimuthal angle of the reaction plane and the brackets denote statistical 
averaging over different events which leaves non-zero only the contribution of the cosine in 
Eq. (0). The reaction plane is determined by averaging over particles produced in a heavy ion 
collision with certain weights designed to optimize the reaction plane resolution Let us 
imagine that the particle with momentum pt in Eq. (|l]) was produced in the same subcollision 
with some other particle, which contributed to determination of the reaction plane angle 
Then the particle pt would be correlated with the reaction plane by these non-flow azimuthal 
correlations. The averaged over events correlations do not disappear at high px, as will be 
shown below. Therefore the elliptic flow observable V2 may be potentially sensitive to the non- 
flow correlations originating from minijet production in the initial conditions. To avoid this 
problem the authors of |Q (see also 0]) introduced a cumulant approach to flow analysis. 
In terms of the saturation model contribution of minijets to higher order cumulants defined in 
^ is suppressed by powers of S±Q'j. ~ Npart- This led the authors of [^, Q to suggest 



that the high-cumulant flow analysis would be relatively free of minijet effects compared to 
standard flow analysis. 

There are certain potential dangers in this conclusion which have to be addressed in an 
explicit calculation. One is that for peripheral collisions S±Ql ~ Npart is not such a large 
number and therefore minijet effects should not be suppressed anymore. From the experimental 



data |21|, ^ we know that elliptic flow is strongest in peripheral collisions making the minijet 
contribution also large. Another potentially dangerous question is whether the parameter S±Ql 
appears in actual calculations with some numerically small prefactor making inverse powers of 
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S±Ql not too small. After all the observed differential f 2 is also not a very large quantity being 
of the order of 10 — 15% at RHIC and may be very sensitive to such corrections. 

To clarify the questions mentioned above we are going to perform an explicit calculation 
of the minijet contribution to the elliptic flow observable. The paper is structured as follows. 
In Sect. 2 we show that the standard determination of flow employing Eq. (|1^) and the flow 
extracted from the 2nd cumulant proposed in |32, 33, 34 1 should yield the same result for V2 
for a trivial choice of weights. This result was confirmed by an explicit flow analysis using two 
different methods at STAR giving almost the same result for V2 • 

The second order cumulant from P^l, 53 is nothing but a two-particle correlation function, 



a contribution to which from the minijet production is calculated in Sect. 3. To construct a 
model of minijet production which incorporates both hard and soft physics we have used 
factorized expression for two-jet production with the unintegrated gluon distributions given by 



the quasi-classical Glauber-Mueller expression H, H, H, gl], m. These gluon distribut 10ns 



are characterized by the saturation scale Qs ^ J^qcd which insures applicability of small 
coupling approaches down to rather small transverse momenta of the produced particles. Our 
model of course does not yield us the exact two-particle production cross section, but gives a 
realistic approximation similar to the one used in [l^, |13| to describe the multiplicity distribu- 
tions at RHIC. An exact calculation of double inclusive minijet production in the quasi-classical 
framework appears to be rather complicated and is left for further research: even the single 
inclusive cross section has not been yet unambiguously theoretically determined, despite the 
extensive efforts |^, |[. 

In Sect. 3 we calculate the differential flow observable resulting from these two-particle 
correlations. Our final result is given in Eq. (0). The obtained differential elliptic flow V2{pt) 
starts increasing as a power of pr for small pr <C Qs (see Eq. (^)) and then saturates to a 
slow logarithmic growth ior px ^ Qs (see Eq. (02)). We also derive the centrality dependence 



of the minijet flow contribution (see Eq. (^)). 

In Sect. 4 after making some simple assumptions about the gluon distributions employed in 
the minijet production cross section we fit the differential elliptic flow STAR data |26| using the 
flow from Eq. (^) with as = 0.3, Qs = 1 GeV^. These values are in agreement with the ones 
used in the saturation- inspired analysis of the multiplicity data in [|1^, [1^. The fit is shown in 
Fig. ^. We can also fit the centrality dependence of V2 with our minijet model by noting that 
the integrated flow should scale as V2{B) ~ 1/ \J Sx_Q1 ~ 1/ Npart (see Fig. 

We discuss the results in Sect. 5 by stating that while a more involved numerical analysis is 
still needed to analyze the emerging RHIC data on elliptic flow we have demonstrated that the 
contribution of minijets to the standard flow analysis is very large, possibly accounting for most 
of high-pT- data. Thus it appears that the standard flow analysis is heavily "contaminated" by 
minijets which prevent direct measurements of the contribution of collective QGP effects to 
elliptic flow. At the same time flow analysis seems to be rather sensitive to details of saturation 
physics and could be used for determination of nuclear saturation scales. 



2 Different Methods of Flow Analysis 



In p2[ Borghini et al proposed a new and interesting approach to flow analysis employing 



higher order cumulants. Let us outline some important features of the approach presented in 
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52| for the lowest order cumulant. The second order cumulant is defined for two particles with 
azimuthal angles 0i and 02 as [B^, I 



,2i(</,i-02) 



,2i(«ii-(/.2) 



(2) 



The average e '^^ y vanishes since the angle here is measured in the laboratory and heavy ion 
collisions are azimuthally symmetric after averaging over many events. For fiow correlations this 
means that the angle in e^*'^^ ^ is not measured with respect to reaction plane, but with respect 
to some fixed direction in the detector. Assuming that the particles 1 and 2 are correlated with 
each other only through the fiow correlations with the reaction plane the authors of |Q wrote 



,2i(0i-02) 



32j(<^l-02) 



^2 



(3) 



where the definition of elliptic fiow from Eq. (|Tp was employed. A new method for measuring 
elliptic flow was proposed in |^ using the two particle (and higher order) cumulants of Eq. (H). 
If we fix the transverse momentum of particle 1 to be px and average over all momenta of the 
particle 2 over many events then as one can see from Eq. (^) 

z^'^^^^P-^^-^^^ ) = v^ipr) {V2) (4) 

with (^2) the elliptic fiow variable averaged over all pr- At the same time if we do not impose 



= {V2). (5) 
noting that only cosine components of the exponents survive the aver- 



any restrictions on the transverse momenta of both particles we get |^ 
From Eqs. and 

aging we obtain the following expression for differential elliptic fiow |^ 

(cos(2(0,(pr)-02))) 

V2{Pt) = , 

y'(cos(2(0i-02))) 

Let us define the event-averaged two-particle multiplicity distribution function 

dN 



(6) 



d?ki dyi d'^k2 dy2 



(7) 



where the transverse momenta of the particles are ki and ^2; while yi and y2 are their rapidities 
and we average over all events with the impact parameter B_ between the two nuclei. This 
distribution can be written as a sum of uncorrelated and correlated terms 



P{k^,k2,B) 



dN 



dN 



+ 



dNr. 



d'^kidyi d'^k2dy2 d'^kidyid'^k2dy2' 



where we suppressed the impact parameter dependence. The correlated term in Eq. 
usually much smaller than the uncorrelated one as it is suppressed by a power of Q'^Sj_ 
Using Eq. (J^) we rewrite Eq. (||) as 



(8) 



IS 



N 



part- 



V2{k 



J d'^k2dy2d(f)idyiP{ki,k2,R) cos(2((/)i - ^2)) 
/ d'^k2dy2d(f)idyiP{ki,k2,R) 
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^ / J d^hdy,d'k2dy,Pihk2,B) 

[j d^hdy,d^k2dy2P{h,k2,B)cos{2{(f),-(j)2))) ' 

where we have relabeled the transverse momenta of the particles to be ki and k2. Eq. (P) 
gives us a way of calculating elliptic flow from two-particle correlation functions. Higher order 



cumulants would yield us ways of calculating flow from higher order particle correlations |32 



though one can not go to arbitrary high order cumulants due to lack of statistics there. An 



analysis of RHIC data has been performed at STAR using both conventional and cumulant 
approaches explicitly demonstrating that the flow obtained by the second cumulant technique 
of Eq. (^) is consistent with the flow extracted using conventional techniques. That is the 
approach of Eq. is equivalent to the conventional flow analysis. To see this let us first note 
that in the actual standard flow analysis one has to take into account the resolution in the 



event plane determination ||3^ . This would modify Eq. (|l]) giving |3T|, |2T 



V2{Pt) = ^^^^ 

/2(cos(2(vl/?j-vI;^)) 



where \E'^ and \E'5? event plane angles determined in two different sub-events labeled 

a and b with the particle multiplicity A^/2 in each of them while is the event plane angle 
determined in the full event with multiplicity N [^. Eq. (p!0| ) is quite similar to Eq. (|^). To 
show that the two equations are almost equivalent let us go back to exponential notation and 
rewrite the numerator of Eq. (^) as 

(cos(2(0i(pt) -02))) = //e'^(MPT)-MkT))\ \ (n) 

\^ ' "■^'■'^Tl^VT / events 

where we first average over all particles (with various k^) in a given event other than the chosen 
particle with momentum pxi just like in the standard flow analysis |2^, and then average over 
all events. On the other hand, the azimuthal angle of the reaction plane ^ r in Eq. ( |T0| ) can be 
defined by 

where Qn is a real number. If one neglects correlations between the particles and multiplies 
Eq. (|12|) by its complex conjugate avergaing over events with the same total multiplicity one 
easily gets < > = 1/A^ |22|. Therefore at the leading order in A^ we have (for large A^) 
Qn ~ I/Vn. In Eq. (|1^) we for the moment forget about the transverse momentum cutoffs 
imposed in event plane determination and subtleties related to different choices of weights and 
detector imperfections. Substituting Eq. (|12|) into the numerator of Eq. (|T0|) we obtain 



Q2i{<f>PT-<t>kT) \ ) (13) 

{Qn) V / all kr^PT/ events 

where we assumed that for a fixed impact parameter collisions the particle multiplicity A^ (and, 
therefore, Qn) is independent of two-particle correlations. Eq. (0) is identical to Eq. (p!l| ) up 
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to a factor of 1/ {Qn)- Similarly one can show that the denominators in Eqs. (|T0|) and are 
equal to each other up to a factor of w2/ 



events 



r)2 \ / all dif ferent kTi,kT2 , 

I events 



-— ^ (cos(2(0i - 02))) . (14) 
Recalling that for large multiplicities {Qn) ~ we see that at the leading order in A^: 



Q\l^l \f2 (Qn) ~ (Qn/2J I {Qn) ~ 1- Extra factors in the numerator and denominator 

of Eq. ( p!0D cancel reducing it to Eq. (H). Thus we showed that Eqs. ( p!0D and (|^) are identical 
in the limit of large multiplicity N. We have proven that the standard definition of flow from 
Eq. (0) is equivalent to the definition introduced in 34 1 given here by Eq. (^). Therefore in 



order to calculate the differential elliptic flow V2{pt) we need only to calculate the two-particle 
correlation function P{ki, k2,B_) and substitute it into Eq. @. 



3 Calculation of V2 

To calculate the two-particle correlation coefficient in Eq. (|^) one needs to know single and 
double particle multiplicity distributions. While the quasi-classical single particle production 
mechanisms in heavy ion collisions have been extensively studied @, H, the double inclusive 
particle production has not been calculated yet. One needs to calculate production of a pair 
of particles with comparable transverse momenta lA;^! ~ 1^2! ~ Qs and rapidities yi ~ y2 in 
the framework of McLerran-Venugopalan model 0, that is resumming all powers of Ql/k^ 
IQ, ^ ^. Unlike the single gluon production this process involves an extra gluon and can not 



be described by the classical field methods of ^ ||. At the lowest order in Ql/k^ the two 
gluon production amplitude is equivalent to the real part of NLO BFKL kernel ^ and is 
rather sophisticated |Q. Going beyond leading order in Ql/k^ appears to be tremendously 
complicated and we will not address this problem here. Instead we are going to construct a 
model of correlated two-particle production employing quasi-classical gluon distributions from 
1^, ^ in the production formula inspired by coUinear factorization ^ , similar to how it was 



done in |]T2|, [T^ in describing the multiplicity data at RHIC. We would also assume that in any 
given event either yi ^ y2 or yi <^ ^2- Of course this assumption does not hold in actual flow 
analyses 23, 21] and we are making it just to simplify the calculations. While, as discussed 



above, a more detailed calculation would still be required to obtain an exact expression for the 
correlation coefficient in Eq. (^, we believe that it would only introduce numerical corrections 
to our approach leaving qualitative results the same. 

Consider inclusive production of one or two gluons in scattering of two quarks at high energy 
as shown in Fig. where the blobs denote effective Lipatov vertices. It is convenient to use 
the Sudakov decomposition of gluons' four-momenta: 

Qi = ttipi + PiP2 + q., i = 1,2,3. (15) 
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Figure 1: (a) One- and (b) two-gluon production amplitudes. Thick dots denote Lipatov 
vertices. 



In Fig. |I](a) the dominant contribution stems from the multi-regge kinematical region where 
ai ^ a2 and (32 ^ Pi- Using Ward identities the t-channel gluon propagators can be written 

as 



s ai(3i ' 



(16) 



where s = 2(pi -^2) is center-of-mass energy squared. With this representation of the t-channel 
gluon propagator effective Lipatov's vertex |^ for s-channel gluon production (black blobs in 
the Fig. P reduces to the usual three-gluon vertex F^^ and can be easily calculated: 



(17) 



In our calculation we will need the square of Lipatov's vertex given by 

(L^r = a,P2s£ql. (18) 
Using Eq. (|l^), Eq. {\17\} , Eq. ([T8|) and the four- dimensional momentum element decomposition 



cPq. dui dPi 



we evaluate the single and double distribution of the produced particles shown in Fig. |1| 



dN 



2 Cf 



d^ki dyi 7r2 Sj_ 



(fqi 



d^kx dyi (Pk2 dy2 n'^ S_l kl kl 
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d'^qi 



ql ih + k^-q^) 



(19) 

(20) 
(21) 



where we changed variables to ki = q^ — q^-, k^ = q^ — q^ (the momenta of produced gluons). 
Throughout the paper we assume for simplicity that the nucleus has a cylindrical shape with 
radius R and "height" 2R and the axis of the cylinder coincides with the collision axis. S± = 
S±{B.) is the nucleus overlap transverse area which depends on nuclei relative impact parameter 
B_. We assume that colliding nuclei are identical having atomic number A each. 

In fact colliding quarks in Fig. are not free but confined inside nucleons. If momenta of 
produced gluons are higher than all scales characterizing the nucleus (including Q^) we can use 
fc^-factorization to write the distributions from Eq. ( ^OD and Eq. (^I]) as convolutions of hard 
amplitudes with parton distributions in nuclei which are predominantly gluonic in the high 



energy region s ^ t. Employing the relation [38 



xGA{x,q^) = AxG{x,q^) = A ln(gV/i2), (22) 

— _ ^ _ 

with yU some non-perturbative cutoff the one and two-particle distributions can be written as 



45 



dN 2 as 1 f ,2 dxGA dxG 



^ 51-?^ 377^^^ (23) 



d^kidyi Cf S± k^ J dq"^ d{ki — q_^ 

dNcorr Nc-ai 1 r o dxGA dxG 



d^ki dyi d^k2 dy2 vr2 Gp S± ^f^^ J ^ dql d{ki + k2 - q^f ' ^'^^^ 

As was argued in the Introduction, if the number of partons in each nuclei is large enough 
then their mutual interactions must be taken into account. This may happen either due to 
increase of parton density in each nucleon by subsequent emission of gluons in course of quantum 
evolution as collision energy increases p|, |lT| or due to enhancement of nucleon parton density 
in a large nucleus by the atomic number A [H, ^ The latter effect leads to creation of the 
non-Abelian Weiszacker- Williams (WW) field A^^{z) of the nucleus which is a strong classical 
gluon field already at moderate energies 0, |[. The WW field gives rise to the unintegrated 
gluon distribution given by 



dxGA{x,q^) 2 



rfg2 (27r; 



2 



rfze-'^l I /6TrU»"»'(0)4"'"'(l)) 



= fi'^^-''-^{l-e-i-"''), (25) 

7r(27r)"^ J ^ ^ 

where b is the gluon's impact parameter (which we can trivially integrate over in a cylindrical 
nucleus) and 

QlU) = ^^^pa;G(x,l/z^)T(6), (26) 

c 

with p = A/[27rR^] the atomic number density in the cylindrical nucleus with atomic number 
A, and T{b) the nucleus profile function equal to 2R for the cylindrical nucleus considered here. 
This provides the initial condition to the nonlinear quantum evolution of the gluon distribution 
with energy in the high parton density region 0. is a scale at which nonlinear nature 
of the gluon field becomes evident. We suggest using the classical expression Eq. (|25|) as an 
approximation to the exact gluon field of the nucleus. This is a justified approximation as long 
as a^lns < 1, i.e., when corrections due to quantum evolution are small. 
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It is usually assumed that fcj-factorization holds in high parton density regime as well as 
in the linear one |^. Phenomenological models for heavy ion collisions which employ this 
assumption together with nonlinear evolution for gluon distributions proved to be successful 
in describing experimental data at SPS and RHIC ||13|- This implies that the use of k^- 



factorization is a quite good approximation for the SPS and RHIC kinematical regions even 
for kx < Qs- Therefore substituting Eq. ( pSf ) into Eq. (p3D and into Eq. (^If ) we obtain the 
following expressions for the single and double gluon distributions 



dN 



CfSi_ AKi dz 



d^kidyi a ski 



Mhz) (i 



— e 



dNr. 



d\ 



(1- 



(27) 



(28) 



d^ki dyi d'^k2 dy2 ki 

where the saturation scales in both nuclei are the same since the nuclei are cylindrical and 
identical. Generalization of Eqs. (^) and ( pSf ) to a spherical nucleus is straightforward. 

We have to point out again that the single particle distribution is, in principle, known better 
than displayed here in Eq. (^) |]^. For instance the distribution in Eq. ( P7| ) is not infrared safe, 
while the correct distribution derived in [^j has no infrared divergences. However, as we need 
both single and double particle distributions to obtain elliptic flow using Eq. (^) we should 
calculate both of them in the framework of the same model. Since the exact calculation of 
the double gluon distribution does not seem feasible at the moment we have to calculate both 
distributions in the same /cT-factorization approach inserting i^-factors Ki and K2 to correct 
the normalization of the approximation to include higher order corrections |T^, which we have 
done in Eqs. (P?]) and (PH]). The value of the i^-factors will be fixed later. We will determine 
Ki by comparing particle multiplicity per unit rapidity [dN/ dy) resulting from Eq. (|27|) to the 



total multiplicity observed at RHIC |T^. To fix K2 we consider two-particle production cross 
section, which is proportional to the two-particle multiplicity distribution function P{ki, k2,B_) 
from Eq. (^. In the limit of large transverse momentum fci ~ ~ Pt the second term in Eq. (^ 
given by Eq. (^) falls off at most as l/p^ and has a collinear singularity fci -|- ^2 = 0- Therefore 
it dominates over the first term in Eq. (^) given by Eq. ( PTj ) squared, which gives a fall 
off of the cross section. For large ki = —k2 Eq. (|^) and, therefore, Eq. (^) should match onto 
the corresponding collinear factorization expression for back-to-back jets [^, ^ . 
To fix the normalization we expand the term in the parentheses of Eq. (|28| ) to the lowest order 
and integrate over one of the transverse momenta obtaining 



dar 



dp^ dyi dy2 



9tt al 



[xG{x,pl)f 



(29) 



41, 46, 47 . We put A^c = 3 



47 are rather successful 



in agreement with collinear factorization result at mid-rapidity [ 
explicitly in Eq. (|29|) . Collinear factorization models 0, 
in describing the high-p^- particle spectra in hadronic and heavy ion collisions when putting 
= 2 to correct the lowest order perturbative expression for next-to-leading order effects. 
For our model to be in agreement with these high-p^ data we will have to also put K2 = 2 
when trying to describe the flow data in the next section. We have to note that since the gluon 
distributions in Eq. ( |25D do not include DGLAP evolution in them and Eq. ( PBD does not 



have jet quenching effects |27, 2S, 47| in it our model can not be applied to describe the high-p^ 
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spectra and a complete DGLAP-evolved gluon distributions along with jet quenching should 
be used in a more detailed numerical treatment of the problem as was done for particle spectra 
in [I 



41, 46 



Substituting Eqs. (^) and ( pS]) into Eq. (||) we are ready to calculate f 2(^75 B_)- To calculate 
the integral in the numerator of the first line of Eq. we note that only the correlated part 
of P{ki, k2,B_) (the second term in Eq. (P)) contributes there. Let us denote the angle between 
two- vectors z and ki by 9 and the angle between ^2 ki by 02- The third angle defining 
direction of k^ is free, so we can just integrate over it obtaining a factor of 27r. Plugging Eq. (|28|) 
into Eq. (|^) and using the latter in Eq. we get 

(fk2d(f)iP{ki,k2,B)cos{2{(f)i-(j)2))= [ d^/c2 (i0i , cos(2(0i - 02)) 



d'^ki dyi d'^k2 dy2 

N,CfS^K2 f dz^^gd^ (^_^^,^Ql/,Y r d^2e''''-^°<f--'^e-^^'^^-^' cos'. 



^ 71^ J z^ kl 
NcCfS±AK2 fdzdklf^ ^z^n^ux'^ 







A;2 7r3 



/5^(l-e-^«^/f J2(M^2(M 



N.C,S.AK2 rdzf,_,-.^QV^)^ j^^k^.y (30) 



k^ TT-^ JO z"^ 

To calculate the integral in the denominator of the first line of Eq. (^ we first note that the 
final state multiplicity per unit rapidity in heavy ion collisions was calculated in |^ (see also 
0) and reads 

_ Cf Ql 
dy a.2 7r2 

with c the gluon liberation coefficient which was calculated in to be c = 2 In 2. Following the 
same steps as in derivation of Eq. ( pOD we find contributions to the integral in the denominator 
of the first line of Eq. (|]) of uncorrected two-particle distribution 

,2, , dN dN dN dN 
d k2 dip 



c ■ 



d'^kidyid'^k2dy2 dy2 kidkidy 
CISIQIAK, foo dz ^^^^ . „-.^q?m2 



a'^kl vr^ Jo z^ 
and of correlated two-particle distribution 



Mkiz) (1 - e-^'«'/^) (32) 



d'^k dip dNcorr 

d^ki dyi d'^k2 dy2 

As one can obviously see the integral in Eq. (^) is enhanced by factor of S^Ql/a1 as compared 
to the integral in Eq. ([53|). This means that the total number of correlated pairs of particles is 



negligible compared to the number of the uncorrelated pairs given by the square of the total 
multiplicity. Therefore the correlated contribution form Eq. (|33D can be neglected compared to 
the uncorrelated contribution in Eq. (p3). 
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Figure 2: Nuclear collision in the transverse plane. 



To estimate the integrals in the second hne of Eq. @ we need to integrate Eq. (^) over 
momenta ki 

J d^h d^k2 Piku k2, B) cos(2(0i - ^2)) 

where we assumed that of Eq. (p6D is approximately independent of transverse coordinate 
z_ neglecting the logarithm of Eq. (|2^) . To complete the calculation of the integrals in Eq. (P) 
we note that ^ 



I d'k^d'k,P{k„k„B) = ^^= (c 
J dyi dy2 \ 



2 1- (35) 

Inserting Eqs. (|30D , (|32D , ( ^If ) and (^) into Eq. (P) and noting that the integration over 
rapidities is trivial and cancels out between different integrals we obtain 

We have changed the transverse momentum back to pr to comply with conventional notation. 
Eq. ( PB| ) gives the minijet contribution to the differential elliptic flow for collision of two identical 
nuclei at given impact parameter B_. V2{pt) measured in experiments is averaged over all impact 



parameters 22 



I d^Bv2{pT,B) {dN/d^pTdy) 
''^^^^ = Jd^BidN/d^prdy) " ^'^^ 

In our model each nucleus is cylindrical, so that the dependence on the impact parameter B 
comes only from the nuclear overlap area S±{B). The overlap area is 

= R^{(3 - sin (3) = (2 arccos(5/2i?) - sin(2 arccos(5/2i?)), (38) 

where (3 is the opening angle in the transverse plane as shown in Fig. ^ From Eq. ( pTf ) we 
see that dN/d^pxdy ~ and from Eq. ( p6l) one concludes that V2{pt,B_) ~ \JS± {B_) . 

Therefore after using Eqs. (|27| ) and ( p6| ) in Eq. (pTf ) the impact parameter averaging reduces 
to 



/o^ 6/5 5^2 arccos(5/2i?) - sin(2arccos(5/2i?) .99 ,„ ^^^^ 



R Jo dB B (2 arccos(5/2i?) - sin(2 arccos(5/2i?)) R 
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where = vri?^ is the cross sectional area of the nucleus. The final result for V2{pt) reads 



/4 



2 



Eq. (^) is our main result for differential elliptic flow. Let us flrst study the qualitative features 
of the flow from Eq. (^OD before using it to flt RHIC data in Sect. 4. 

In the small momentum region, pt ^ Qs, we can neglect the exponents in the numerator 
and denominator of Eq. Expanding the Bessel functions for small pxz and cutting off the 
integrals over z hj 1/pt from above and 1/Qs from below we end up with 

We see that minijet f 2 is an increasing function of transverse momentum, which is in qualitative 
agreement with the data. 

The asymptotic behavior of V2 at high transverse momenta of the particles, pt ^ Qs, can 
be found by expanding the exponents in the numerator and denominator of Eq. (|40|). Here it 
is essential to keep logarithms arising from Eqs. (^Gf) and (^). Simple integration yields 

(vr^ N K \ ^^'^ 



Therefore we conclude that while the differential elliptic flow given by Eq. (40) increases for 
small transverse momenta, at pt ~ Qs it turns over and its growth slows down to just a 
logarithmic increase. As the transverse momentum increases the DGLAP logarithms should 
become important modifying Eq. (|42|). The logarithm in Eq. (^) might also be an artifact of 
our approach and it might disappear if a more detailed analysis is carried out with inclusion of 
evolution effects in the gluon distribution functions. At the moment we can make an observation 
that the contribution from minijets to differential elliptic flow is in qualitative agreement with 
the data PB |. 

To understand the centrality dependence of the elliptic flow coming from minijets we use 



the deflnition 22 



^r^^ ^ / d'^PT dy V2 {pt, R) {dN / (Ppr dy) 
""^^ ' J d^pTdy{dN/d^pTdy) ^ ^ 



to obtain using Eqs. ( ^Tf ) and (|3T1) 



/ 21n27r^iy, 

V2{B) = a, o /oN^2 ' (44) 



yC^CFS^{B)Ql^ 

where Sj_{B) is now the 5-dependent nuclear overlap area. Since Sj_{B) Ql ~ Npart we may 
therefore conclude that minijets flow scales as 

V2{B) ~ -^=. (45) 

Elliptic flow from Eq. (^5]) is smaller for central collisions with large Npart and it increases 
for peripheral collisions with decreasing Npart, again in qualitative agreement with the data 
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4 Our model versus experimental data 



There are three interesting kinematical regions in high energy heavy ion colhsions corresponding 
to various values of transverse momenta pt'- 

(i) Pt^Qs where multiple rescatterings of partons dominate. In this region transverse mo- 



mentum dependence of V2 is given by Eq. (^Tj). Owing to the nonlinear interactions nuclear 



structure functions in region (i) are functions of only one variable Q/Qs{x) instead of two 



Q and x. This is often referred to as geometric scaling [50 



(ii) Pt ^ Qs where the DGLAP evolution starts to be important. There the leading 
behavior of V2 is given by Eq. ( ^2] ) in which DGLAP equation may introduce logarithmic 
corrections, possibly changing the power of ln(p7^//i). The scaling behavior is broken 
down. 

(iii) At high energies in the course of quantum evolution the saturation scale Qs becomes much 
larger than any soft QCD scale A. It was argued both analytically and numerically 
jSOl , |To| , |52| that the geometrical scaling holds in a much wider kinematical region than 
(i). In fact the additional scaling region is Qs^Pt^QH ^- Here even though the scatter- 
ing amplitude is still far from saturation, the nonlinear interactions produce remarkable 



scaling phenomenon. It was argued in [|T2[ that the saturation scale in RHIC kinematical 
region is of the order of 1 ^ 1.4 GeV. Since A ~ 0.1 ^ 0.2 GeV, we expect the scaling to 
be important at p-r < 5 ^ 10 GeV, i.e., throughout most of the kinematic region of the 



differential elliptic flow data ||26 



To compare predictions of our model with data collected by the STAR collaboration at 
RHIC ||2^ we should evaluate f 2 given by Eq. (^) and Eq. (|^) at < 5 GeV. This means 



calculation in the regions (i) and (iii) where geometric scaling works and the only relevant 
dimensional scale is Qg. Quantum evolution changes the quasi-classical result of Eq. (|26|) so 
that the saturation scale acquires energy dependence. Its analytical expression can be calculated 
from the nonlinear evolution equation 0. There is a simple way to introduce Qs into our model 
preserving geometric scaling developed by quantum evolution, where Qs will play a role of a 
phenomenological parameter which is used to fit the data. Namely, by definition the saturation 
scale is a scale at which the expression in the Glauber exponent of Eq. equals one 



|^^Q^U)l.^^4/gi = ^ln| = l, (46) 

where we have introduced in the quasi-classical (Glauber) approximation 

o Anal A , . 

Qlo = — (47) 

which is an unknown constant for the case of a fully evolved distribution. Qso can be easily 
expressed in terms of Qs using Eq. ( ^BD 

Qso = (48) 

2(1 
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Substituting Eq. (|i8|) into Eq. (50) means the following change in the power of the Glauber 
exponents 

Using Eq. (|49|) in Eq. (^0]) we can write for the integrals involved 

^ Uvrz) (l - e-^^«^(^)/f =vir% 1 - e-(^^«^/^-^) (50) 



where we defined ^ = pt^, A = 2/i and n = 0or2. Geometric scaling implies that the 
distribution functions depend only on one parameter pr/Qs- To eliminate the scaling violating 
terms in Eq. ( |5(]| ) we note that the average value <z>~ 1/ <pr>~ and rewrite the 

logarithm as 

In^^ln^ (51) 

where A would be an independent parameter of our fit restricted by reasonable possible values 
of the non-perturbative scale. The integral of Eq. (^) becomes 



Pt 



oo 



pJ„(0 [l-e"(«^«^/^^-)TiiS) I . (52) 



This form of the integrals will be used in Eq. (^D|) to numerically estimate V2{pt)- 

To obtain numerical value of V2{pt) we have to estimate the normalization coefficient Ki. 
To this aim note that the fraction of the total particle multiplicity per unit of rapidity due to 
soft saturation physics in heavy ion collisions at center-of-mass energy ^/s = 130 GeV is given 

by m 

dNsat ^ ^ 3 A^part .p-oN 

dy ^ ^-^^ 2 ' ^^^^ 

where 3/2 is the conversion factor from charged particles to all particles multiplicity. On the 
other hand, we can calculate the total particle multiplicity by integrating Eq. ([27| ) over all 
momenta ki which gives 

iJi^j,^^^n2SiC,Ql ingj. (54) 

ay TT^ as 2fi 

Equating Eq. ( |53D and Eq. ( |MD for the head-on collisions {B = 0) we obtain 

^ 3 7r^a,1.02Arp,rt(g = 0) 
' 4\n2SiCFQlHQJAy ^ ' 



where A'part(5 = 0) = 344 |T2|. For gold nuclei with a, = 0.3, Qs = 1.0 GeV and A ^ 0.15 GeV 



the normalization factor is Ki ^ 0.14. 

With the help of Eq. (^) the integration in Eq. (pOf ) has been done numerically for A ^ 
0.15 GeV. We put K2 = 2 in Eq. (|40|) in agreement with coUinear factorization approaches 
1^. Sf = TiR^ with the radius given by Woods-Saxon parameterization R = I.IA^^^ fm. 



The results are shown in Fig. ^. It can be seen that our model describes the high pt data 
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Figure 3: Differential elliptic flow data from STAR |2q/ versus the predictions of our niinijet 
model. Different lines correspond to our predictions with different values of the saturation 
scale: Qs = 1.0 GeV (solid line), Qs = 1.1 GeV (dashed line) and Qs = 0.9 GeV (dash-dotted 
line). We used the following parameters A=0.15 GeV, as=0.3, A=197 (Gold). 

remarkably well. We are led to the conclusion that the high momentum tail of the V2 is 
saturated by the two-particle correlated minijet production. We can also explain the turnover 
of z;2 at pt ~ 1 -j- 2 ~ GeV as being due to saturation effects in the wave functions of the 
colliding nuclei at low pt- Another important observation is that two-particle correlations in 
the initial state give significant contribution at low pt as one can see in Fig. Note that the 
p^ increase of f 2 at small momenta may be an artifact of our model which neglects the change 
of the anomalous dimension 7 of the gluon structure function. It is well known that correct 
anomalous dimension at the boundary of the saturation region is 7 ~ 1/2 |T^. It smoothly 
changes from in the unitarity limit to 1 in the Bjorken limit |jlO|. More accurate analysis 



is needed to understand the low pr behavior of V2 due to two-particle correlations in minijet 
production. 

Applicability of our model is restricted to the mid-rapidity kinematical region. However, 
since the total particle multiplicity is dominated by the particle production at mid-rapidity 
the averaging over all impact parameters in Eq. (|37D gives a reasonable approximation. The 
situation is different if we want to apply our model to describe the centrality dependence of 
V2. While we expect it to be in agreement with the data for the most central events, it is not 
legitimate to use it for peripheral collisions. 

Eqs. and (^) imply 



'Np^rt V 2 In 2 Cf 6^ Qj Np^n J 



where we used c = 2 In 2 . In Fig. ^ we show the result of numerical calculations compared 
to STAR data |26|. The data is represented by black dots while our fit is given by empty 



crosses. We have used Eq. ( p6D with the coefficient K2 = 2. To find the average A^part in a 
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Figure 4: Centrality dependence of elliptic How given by the STAR data (black dots) and our 
fit (empty crosses). We used A=0.15 GeV, a,=0.3, A=197, Qs=l GeV and Np.^,tiB = 0) = 
344 111/. 



given centrality bin we have used the model of Kharzeev and Nardi ||T2 
in our fit correspond to the widths of centrality bins 



Horizontal error bars 



Vertical error bars account for our use 
of cylindrical nuclear shape instead of spherical and for the uncertainty in A^part coming from 
the model of |T^. As expected our model provides a reasonable description of the data for 
central events while the fit does not work that well for peripheral collisions. The saturation 
scale decreases toward the large impact parameters where at some point it becomes of the order 
of Aqcd and the small coupling approach breaks down. Thus, as was mentioned above, our 
approach is not applicable to very peripheral collisions and can not be expected to work for 
them. We conclude that our model gives a reasonable description of V2 centrality dependence 
for the values of centrality nch/'^max > 0.4. 

The maximum of the centrality distribution of elliptic flow is probably determined by the 
onset of non-perturbative effects. The flow from Eq. (44) increases with increasing B, until 
at some point the saturation scale Qs{B) becomes of the order of Aqcd and non-perturbative 
effects take over keeping V2{B) approximately constant for peripheral collisions. (Strictly speak- 
ing Eq. (^^ was derived for cylindrical nuclei with saturation scales independent of the impact 
parameter. However in real life Qs is of course a function of B [^, |T2|, |13| .) Let us denote the im- 
pact parameter at which the non-perturbative effects take over by Bq, so that Qs{Bo) ^ Aqcd- 
Since Qs is an increasing function of the center of mass energy s [y, the impact parameter 
Bq must also be an increasing function of s. That is as energy increases more and more collisions 
become perturbative with the non-perturbative physics being responsible only for increasingly 
more peripheral collisions. The maximum of the flow centrality distribution could be estimated 
from Eq. ( ^ to be V2{Bq) ~ 1/ ^S±{Bo) Aqqjj. With the increase of energy Bq increases, 
S±{Bq) decreases and, therefore, V2{Bq) increases, which qualitatively agrees with RHIC and 



SPS data 26, 19, EO, 21, 22, 23 



16 



5 Conclusions 



In this paper we have calculated the contribution of pairwise azimuthal correlations in mini- 
jet production to elliptic flow variable. The resulting differential elliptic flow V2{pt) given by 
Eq. (|40|) is in good qualitative (see Eqs. ( PD and (|42D ) and quantitative (see Fig. ^ agree- 
ment with the emerging RHIC data pT| , p2| , |23| , |26| . The centrality dependence of the minijet 
contribution to elliptic flow V2{B) (see Eqs. (^^ and (|45|) ) successfully describes the data for 
sufficiently central collisions as shown in Fig. H. The maximum V2{B) appears to increase with 



energy also in agreement with the data |26 



There are several important questions which still have to be addressed in the future as more 
flow data is produced at RHIC. One question concerns the value of the directed flow Vi. It may 
seem from the above discussion that using two-particle correlations one might write, similar to 
Eq. (I), 

/ N (cOs(0i(pt) -02)) 

MPt) = — , =^ (57) 

.\/(cOs(0i - 02)) 

and using the correlations of Eq. (^) get a non-zero directed flow. However, let us remind 
the reader that in the analysis of directed flow the signs of the weights used in determination 
of the reaction plane are reversed in backward hemisphere with respect to the forward one 
Tsj |3T| . That is, in the numerator of Eq. (^7|), the contributions from y2 > and ?/2 < 



come in with different signs (we denote the central rapidity of the event by y = 0). Therefore, 
since in our boost-invariant model the two contributions are identical, the final result for vi 
at mid-rapidity is zero in agreement with the SPS data If we impose some constraint on 
the rapidity interval of the analyzed particles making it asymmetric with respect to y = 0, 
e.g. analyzing only particles with y > 0, the differential directed flow from minijets would also 
become non-zero as observed experimentally ||19|| . 

In our analysis we have neglected contributions from correlated production of three and 
more particles in a single subcoUision. The corresponding diagrams may also contribute to the 
two-particle correlation function in Eq. (H). Of course they are suppressed by extra powers of the 
strong coupling constant compared to Eq. (ESf), though these extra powers are compensated 



by the large logarithms arising due to phase space integration P8| , In heavy ion collisions 



at RHIC the saturation scale is of the order of ~ 1 GeV |13], ^ and the corresponding 
value of the coupling constant is as{Qs) ~ 0.3, which is not too small and will also be enhanced 
by logarithms of energy and transverse momentum. Thus to have a precise description of the 
minijet contribution to elliptic flow one has to resum all these multiple emission diagrams, 
which is equivalent to including the effects of quantum evolution in the double inclusive gluon 
distribution of Eq. (p8|). Since at the moment there is even no rigorously derived expression for 
double gluon production in AA in the quasi-classical approximation, including evolution effects 
in it seems like an important but not immediate goal. 

To improve the quality of the minijet predictions one has to also relax the yi ^ y2 (or 
yi ^ 1/2) condition which we imposed to simplify the calculation. This condition led to absence 
of rapidity correlations in the two-particle distribution given by our Eq. (0). Also the azimuthal 
correlations given by Eq. ( |2^ ) are only "back-to-back" , i.e., the produced particles are correlated 
mostly in the opposite directions in the transverse plane. This behavior seems to contradict 
recent results reported by PHENIX E^, |53| where the two-particle correlation function C ( A0 = 
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01 "~ 02) has maxima at both A0 = tt (back-to-back) and A0 = 0. Relaxing the yi ^ 1/2 (or 
yi ^ 2/2) condition would greatly complicate the particle production calculations though 
would give predictions for elliptic flow in a much more realistic kinematics. The appropriate 
calculations were performed for production of a pair of particles with i/i ~ 1/2 by Leonidov 
and Ostrovsky in The results of |Q indicate that the two produced particles are indeed 
somewhat correlated in rapidity. The analysis of azimuthal distributions performed in |44] also 
demonstrates enhancement of both A0 = vr and A0 = correlations in the obtained double 
particle spectrum, in qualitative agreement with PHENIX data These correlations 

should also contribute to smallness of directed flow Vi. A detailed analysis of elliptic flow 
employing the full calculation of |Q will be done elsewhere |5^. In a complete analysis one 



has to put in realistic geometrical shapes for the nuclei as well. Of course the effects of different 
cuts used to analyze the data and various choices of weights have to be incorporated in the 
precise analysis too. While these modiflcations are likely to introduce some numerical changes 
in the flts to the elliptic flow data presented above, we believe that the qualitative conclusions 
of importance of the minijet production to flow variables would remain unchanged. 

We therefore conclude that minijets give a large contribution to the elliptic flow extracted 
using current flow analysis methods, probably accounting for most of the flow at large Pt (see 
Fig. Differential elliptic flow appears to be sensitive to saturation physics of the early stages 
of the collisions. To study QGP it would be very useful to invent a method of flow analysis 
that would be insensitive to minijet contribution and would only measure the collective effects 
due to elliptic flow [0, |3^. We have to note, however, that if one tries to calculate the double 
inclusive minijet production cross section in the saturation approach p|, ^, |^ many of the 
diagrams that would contribute would reduce to independent production of two gluons which 
would then exchange a gluon with each other. This 2—^2 process would constitute the flrst 
step in the onset of thermalization [0. Therefore it appears that correlated production of pairs 
of minijets may be intrinsically related to thermalization. On one hand this implies that the 
minijet flow of Eq. ( ppj ) is partially due to collective effects after all. On the other hand it may 
be impossible to create an observable which would distinguish these early stage thermalization 
effects from the flow of the fully thermalized quark-gluon plasma. 
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